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^propose  an  approximate  Newton  method  for  solving  thf  coupled  nonlinear  system. G(t<,  f)  =  0 
and  jV(u,t)  =  0  where  ti  €  f2",  t  €  G  :  R?  x  Kn  >-+  R*  and  N  :  R"  x  Rm  Rm.  'The  method 
involves  applying  the  basic  iteration  5  of  a  general  solver  for  the  equation  G(u,  t)  =  0  with  t  fixed. 
It  is  therefore  well-suited  for  problems  for  which  such  a  solver  already  exists  or  can  be  implemented 
more  efficiently  than  a  solver  for  the  coupled  system.  We  deriveconditions  for  S  under  which  the 
method  is  locally  convergent.  Basically,  if  S  is  sufficiently  contractive  for  C,  then  convergence  for 
the  coupled  system  is  guaranteed.  Otherwise,  we  shew  how  to  construct  a  S  from  S  for  which 
convergence  is  assured.  These  results  are  applied  to  continuation  methods  where  N  represents  a 
pseudo-arclength  condition.  We  show  that  under  certain  conditions  the  algorithm  converges  if  S 
is  convergent  for  G. 
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1.  Introduction 

In  this  paper,  we  are  concerned  with  computational  algorithms  for  solving  coupled  nonlinear 
systems  of  the  form: 


where  z  £  (u,t),u  €  /P*,t  €  /P",  G  :  Rn  x  Rm  *-►  /?"  and  N  :  ft"  x  /fm  *-►  Rm.  More  general 
coupled  systems  can  always  be  casted  into  this  form.  We  shall  assume  that  a  solution  z*  exists 
and  that  it  is  regular,  i.e.  the  Jacobian 


Gt 

Nt 


) 


is  nonsingular  at  z*.  f 

Since  z*  is  a  regular  solution  of  C(z)  =  0,  many  conventional  iterative  algorithms  can  be 
applied  to  solve  for  z*.  However,  this  approach  may  fail  to  exploit  certain  structures  which  are 
inherent  in  the  operators  G  and  N  but  which  do  not  exist  in  C.  Such  structures  could  be  symmetry, 
positive  definiteness,  separability,  sparsity  and  bandedness.  Exploiting  these  structures  may  be 
crucial  for  the  overall  efficiency  of  the  computational  algorithm,  especially  for  large  problems. 

In  addition  to  these  general  properties,  one  may  have  special  knowledge  of  G  and  N,  perhaps 
already  implemented  in  easily  available  efficient  solvers,  whereas  such  solvers  may  not  exist  for 
the  coupled  system.  Situations  like  this  occur  quite  often  in  applications  to  continuation  methods, 
optimization  problems  and  coupled  partial  differential  equations.  In  continuation  methods,  G  may 
represent  a  nonlinear  system  in  u  with  dependence  on  some  parameters  t  and  N  may  represent  an 
arclength  condition  constructed  to  follow  the  solution  manifolds.  If  G  represents  a  discretization 
of  a  well  studied  mathematical  model  (e.g.  the  Navier-Stokes  equations  with  t  being  the  Reynold’s 
number),  one  may  have  special  solvers  for  G  (for  fixed  t)  whereas  these  special  techniques  cannot  be 
easily  adapted  to  solve  the  coupled  system  [4].  Similar  situations  occur  in  constrained  optimization 
problems,  where  t  may  represent  the  Lagrange  multipliers  and  N  the  constraints.  In  coupled  partial 
differential  equations,  for  example  those  that  arise  in  semiconductor  modelling  [8],  G  (with  t  fixed) 
may  be  some  standard  differential  operator  for  which  special  efficient  solvers  exist  whereas  no  such 
efficient  solvers  exist  for  the  coupled  system. 

For  the  above  reasons,  in  this  paper  we  shall  consider  a  special  class  of  algorithms  for  solving 
the  coupled  system  which  makes  use  of  a  general  solver  (presumably  efficient)  for  G,  for  fixed  t. 
We  shall  assume  that  this  solver  is  available  in  the  form  of  a  fixed  point  iteration  operator  S,  which 
takes  an  approximate  solution  u,-  and  produce  the  next  iterate  u,+j  =  S(u t).  Since  most  solvers 
for  nonlinear  systems  are  iterative  in  nature,  it  should  be  relatively  easy  to  extract  S  from  them. 

We  emphaize  that  it  may  not  be  straightforward  to  incorporate  such  a  solver  S  in  most 
conventional  methods  for  solving  the  coupled  system.  The  most  obvious  approach  is  to  use  S  in 
conjunction  with  a  block  relaxation  method  in  which  S  is  used  to  solve  for  u  as  a  solution  to  the 
equation  G(u,t)  =  0  with  1  fixed.  However,  such  an  approach  will  most  likely  fail  if  the  Jacobian 
J  of  the  coupled  system  is  not  positive  definite  or  diagonally  dominant  near  the  solution,  which 
is  usually  the  case  if  the  coupling  between  the  two  equations  is  strong.  Another  classical  method 
which  has  much  better  local  convergence  properties  is  Newton’s  method.  However,  a  linear  system 
with  J  must  be  solved  at  each  step.  This  requires  some  approximations  to  G«  and  G t  which  may 
not  be  readily  available  (in  5  or  otherwise).  Moreover,  while  it  is  possible  to  exploit  a  solver  for 
Gm  when  solving  for  the  linear  system  with  J  (5,  12],  it  is  not  obvious  how  to  exploit  the  special 
solver  5  if  it  does  not  use  the  Jacobian  G%  explicitly. 


tTfcroaghout  thi*  paper,  tubeeripU  in  «,  (  ud  *  denote  partial  differentiation. 


Thus  ii  seems  desirable  to  have  a  class  of  algorithms  for  solving  the  coupled  system  which 
can  be  proven  to  be  convergent,  at  least  locally,  under  rather  general  conditions  but  which  can 
also  effectively  exploit  a  special  solver  5  for  G.  Such  an  algorithm  would,  for  example,  allow 
continuation  techniques  to  be  easily  applied  to  an  application  area  in  a  modular  fashion  and  with 
the  efficiency  built  into  special  solvers  specifically  designed  for  the  application.  It  would  also  allow 
constraints  to  be  added  to  a  special  solver  for  a  class  of  unconstrained  optimization  problems 
without  a  sacrifice  in  computational  efficiency. 

In  Section  2,  we  present  such  an  algorithm  which  is  based  on  Newton’s  method  for  solving 
the  coupled  system  C(z)  =  0.  The  basic  idea  is  to  use  S  to  approximately  solve  the  linear  systems 
involving  J  at  each  step  of  Newton’s  method.  If  S  represents  one  step  of  Newton’s  method  for 
the  equation  G(u,t)  =  0  with  t  fixed,  then  the  algorithm  reduces  exactly  to  Newton’s  method 
for  the  coupled  system,  with  a  block  elimination  algorithm  [5,  12)  applied  to  solve  the  systems 
with  J.  If  S  implements  any  other  convergent  method  for  G(u,t)  =  0,  then  the  linear  systems 
for  J  are  solved  only  approximately.  In  this  way,  the  algorithm  can  be  viewed  as  an  inexact 
Newton  method  [7],  except  that  the  size  of  the  residual  is  not  directly  controlled.  In  Section  3, 
we  analyze  the  local  convergence  properties  of  this  algorithm.  Basically,  we  prove  that  if  p(S%)  or 
||5«||  is  sufficiently  smaller  than  1,  then  the  algorithm  is  locally  convergent.  In  other  words,  if  S 
implements  a  reasonably  fast  convergent  method  tor  G,  then  the  algorithm  will  converge  locally 
for  C.  If  p(Sm)  or  ||5a||  is  not  small  enough,  we  show  how  to  construct  a  S  from  S  for  which 
convergence  is  assured.  In  Section  4,  applications  to  arclength  continuation  methods  are  discussed. 
We  show  that  under  rather  mild  conditions  the  algorithm  is  locally  convergent  if  S  is  convergent 
for  G.  Some  concluding  remarks  are  given  in  Section  5. 

S.  Algorithm 

At  each  step  of  Newton’s  method  applied  to  C(z)  =  0,  the  following  linear  system 

(S  2)0  -0 

has  to  be  solved  for  the  changes  [Su,  St)  in  the  Newton  iterates.  In  order  to  exploit  a  solver  for  G, 
(assuming  it  is  available),  the  above  system  is  often  solved  by  the  following 

Block  Elimination  Algorithm:  [12] 

1.  Solve  Gnw  ss  -G  for  w,  where  w  €  IP*. 

2.  Solve  G*v  =  G%  for  v,  where  v  €  IP*  x  ft™. 

3.  Solve  (JV,  -  Nnv)6t  =  -(N  +  JV»  for  ft. 

4.  Compute  Su  =  w  -  vft. 


Note  that  m  +  I  linear  systems  involving  G*  have  to  be  solved.  Now  assume  that  we  have  a 
solver  for  G(u,t)  =  0  in  the  form  of  a  fixed  point  iteration  u  «-  S(u,t)  with  t  fixed.  For  example, 
Newton’s  method  for  G  would  correspond  to 

sN«.Um  m  u  _  G;l(U,f)C(u,t). 

The  idea  in  the  new  algorithm  is  to  use  5  to  approximately  solve  fo  w  and  v  in  Steps  (1)  and  (2) 
in  the  Block  Elimination  Algorithm.  Since  the  vector  w  is  precisely  the  change  in  the  iterate  u  in 
one  step  of  Newton’s  method  applied  to  G(u,l)  »  0,  it  seems  natural  to  approximate  w  by 


P*ge  3 


where  u,  and  t,-  are  the  current  iterates.  Hie  situation  for  approximating  v  is  slightly  more 
complicated  since  it  does  not  directly  correspond  to  an  iteration  based  on  G(u,  t)  =  0.  However, 
note  that  by  differentiating  the  equation  G(u,  t)  =  0  with  respect  to  t  we  obtain 

C|Ut  +  Gt  —  0 

and  thus  at  convergence 

V  =  — U|. 

Since  at  convergence  u  *  S(ti,t),  it  follows  by  differentiation  that 

Uf  B  Smui  +  S( 

at  x*.  Thus  if  5  is  sufficiently  contractive  for  G  (for  example  if  ||5«||  is  sufficiently  small),  then 
it  seems  reasonable  to  approximate  v  by  — S*.  In  particular,  if  5  =  sNtwl<m  then  =  0  and  this 
approximation  is  exact.  If  S  can  easily  be  differentiated  with  respect  to  t  (for  example  if  S  is  linear 
in  I),  then  St  can  be  computed  without  too  much  difficulty.  In  general,  St  can  be  approximated 
by  finite  differencing  5  with  respect  to  t.  We  summarize  the  above  in  the  following  : 

Algorithm  ANM  (Approximate  Newton  Method)  :  Given  an  initial  guess  («o,  to),  iterate  the 
following  steps  until  convergence: 

1.  Compute  w  =  S(u, -  u,-. 

2.  For  y=l,m  compute 

S(ui,U  +  (jtj)  -  S(u„t,) 

=  - 

where  vj  denotes  the  j-th  column  of  v,  tj  denotes  a  small  finite  difference  interval  and  tj 
denotes  the  j-th  unit  vector. 

3.  Solve  the  following  m  by  m  system  for  d: 

(N,(ui,ti)  -  Nm(ui,ti)v)d  =  ~(N(ui,ti)  +  Nu(u{,ti)w) 

4.  Compute  tl+|  =  I,-  +  d. 

5.  Compute  «f+|  ««,•  +  »-  vd. 

Note  that  similar  to  the  Block  Elimination  Algorithm,  m+1  calls  of  5  are  needed  per  iteration. 
Moreover,  for  this  algorithm  to  be  well-defined,  (Nt  -  Nnv)~l  must  exist  at  all  the  iterates  so  that 
the  linear  system  for  d  can  be  solved.  For  S  -  Syewt<m ,  we  shall  show  that  (Nt  -  Nuv)~*  does 
exist  at  x*  if  Gm  is  nonsingular  there.  For  it  follows  from  a  LU- factorization  of  J  with  G%  as  pivot 

that 

det(J)  «  det(Gm)det(Nt  -  N%G~lGt); 

and  since  at  x*  det(J)  /  0  and  v  =  -ui  *  G^lGt,  it  follows  that  det((Nt  ~  N*v))  £  0  at  x*. 
Therefore,  for  a  general,  sufficiently  contractive  solver  S  (so  that  t;  »  -u<),  it  is  reasonable  to 
assume  that  (Nt  -  Nmv)~l  exists  locally  around  the  solution  x*. 

S.  Convergence 

In  this  section,  we  analyze  the  local  convergence  of  Algorithm  ANM.  For  this  purpose,  we 
view  the  algorithm  as  a  fixed  point  iteration  F: 

-  F(*j). 


We  note  that  the  necessary  and  sufficient  condition  for  convergence  is  f(Ft)  <  1  at  the  solution  z* 
and  a  sufficient  condition  is  ||F,||  <  1  in  some  norm.  We  shall  denote  the  first  n  components  of  F 
by  Fx  and  the  last  m  components  by  F*.  Since  all  relevant  quantities  in  this  local  analysis  are  to 
be  evaluated  at  the  solution  r*,  from  now  on  we  shall  drop  all  the  arguments.  We  also  note  that 
at  **,  we  have  w  =  0  and  d  =  0. 

To  simplify  the  analysis,  we  shall  write 


where  t,  with 


v  =  -St  + 1 


ll«ll  =  M). 


represents  the  truncation  error  in  the  finite  difference  approximation  of  v  to  -St. 

In  order  to  evaluate  F„  we  need  to  compute  Fj ,  F,1 ,  Fj  and  Ff.  From  the  definition  of 
Algorithm  ANM,  it  follows  that,  at  2*, 

Fj  =  /  +  wM  -  [vd)m  =  /  +  «,-  vdu , 

F,1  =  »i  -  (vd)t  -  wt  -  vdlt 

F*  =  d„, 

Therefore,  we  need  to  evaluate  wM,wt,dM  and  dt  at  z*.  From  the  definition  of  w  in  Step  (I)  of 
Algorithm  ANM,  it  is  easily  seen  that 

m  Sm  —  /, 

v(  =  S|. 

After  some  manipulations,  it  can  also  be  verified  from  the  definition  of  d  in  Step  (3)  of  Algorithm 
ANM  that,  at  the  solution  z*, 

dn  =  ~(Nt  -  N%v)~lNnS%, 
dt**-I-(Nt-Nnv)-lNM€. 

Combining  these  results  gives  the  following  : 

Lemma  S.1.  At  the  solution  z* , 

F  -  (  P5«  P(  > 

*  ■  V -(N*  -  Nuv)~lNnS%  -(AT,  -  Nnv)~lN*€  )  ’ 

where Ps/  +  v(iVj  -  Nnv)~xN%. 

From  Lemma  3.1  it  follows  directly  that  if  ||5S||  and  ||<||  are  sufficiently  small  in  some  norm, 
then  ||F«||  <  1  and  Algorithm  ANM  converges  locally.  If  St  is  directly  computable,  then  t  =  0.  If 
finite  differencing  is  used  and  if  the  variables  are  scaled  appropriately,  then  <  can  usually  be  chosen 
to  be  of  the  order  of  the  square  root  of  the  machine  precision  (10]  which  in  most  cases  is  much  less 
than  1.  Therefore,  to  simplify  the  analysis,  we  shall  take  t  =  0  from  now  on.  This  assumption 
should  have  a  relatively  minor  effect  on  the  local  convergence. 

With  this  assumption,  we  have 


p  m  (  PS •  °\ 

*  \-[Nt-N.v)-'N%Sm  O)' 


Since  s(F«)  »  ?{PSm),  it  immediately  follows  that 


Theorem  8.1.  Algorithm  ANM  converges  iff p(PSm)  <  1. 

Specifically,  if  S  =  SNrwUm,  we  have  Sa  =  0  and  therefore  the  quadratic  convergence  of 
Newton’s  method  for  the  coupled  system  is  recovered.  On  the  other  hand,  it  would  be  nice  to 
determine  sufficient  conditions  on  the  contractivity  of  S  for  Algorithm  ANM  to  be  convergent. 
First  of  all,  we  have  the  following  general  sufficient  condition: 

Theorem  3.1.  Algorithm  ANM  converges  if  ||Sa||  <  pjj,  in  a oy  vector  induced  norm. 

Proof.  Follows  from  p(PS9)  <  ||PSa||  <  ||/>||  ||S«||. 

I 

Since  in  general  p(Sm)  <  ||Sa||,  it  is  desirable  to  have  less  stringent  conditions  on  p(Sa)  instead. 
Unfortunately,  this  is  possible  only  if  P  and  Sa  belong  to  special  classes  of  matrices. 

Lemma  3.3.  If  S9  is  normsl,  then  p(PSM)  <  ||P||2^(5.).  If  in  addition,  P  is  normal,  then 
f(PSm)  <  p(P)p(S9). 

Proof  Follows  easily  from  the  fact  for  any  matrix  A ,  p(A)  <  ||A||j  with  equality  if  A  is  normal. 

I 

Lemma  3.3.  If  P  and  Sm  are  simultaneously  diagonalizable,  and  the  corresponding  eigenvalues  of 
P  and  Sm  are  *,  and  <fj,  then  p{PSm)  =  maxi<j<9 

These  two  lemmas  give  the  following  sufficient  conditions  on  p(S%)  for  the  convergence  of 
Algorithm  ANM: 

Theorem  33. 

IfS%  is  normal,  then  Algorithm  ANM  converges  if  p(Sm)  <  lfP  and  S«  are  both  normal, 

or  are  simultaneously  diagonalizable,  then  Algorithm  ANM  converges  if  p(Su)  <  ^y. 

We  note  that  in  general  p{Su)  <  1  is  neither  a  necessary  nor  a  sufficient  condition  for  the 
convergence  of  Algorithm  ANM.  In  other  words,  it  could  happen  (and  we  have  carried  out  numerical 
experiments  confirming  it)  that  a  non-contractive  S  (or  G  can  lead  to  a  convergent  Algorithm  ANM. 
This  can  happen,  for  example,  if  P  and  Sa  are  simultaneously  diagonalizable  and  P  has  a  small 
eigenvalue  corresponding  to  a  large  eigenvalue  of  5«  (or  vice  versa),  so  that  the  product  is  smaller 
than  1.  In  this  way,  P  can  be  thought  of  as  a  projection  (an  oblique  one  in  general)  operator.  In 
practice,  however,  it  would  only  be  prudent  to  employ  a  contractive  S  with  p(S9)  <  1. 

Theorems  3.2  and  3.3  give  upper  bounds  on  p(Sm)  and  ||Sa||  for  the  convergence  of  Algorithm 
ANM.  Based  on  the  assumption  that  (JVj  -  Nmv)  and  G%  are  nonsingular  at  all  the  iterates,  it 
follows  that  p(P)  and  ||P||  are  bounded.  Therefore,  the  upper  bounds  for  p(S%)  and  ||Sa||  in 
Theorems  3.2  and  3.3  are  bounded  away  from  zero.  The  size  of  p(P)  and  ||P||  depends  on  both  N 
and  5  and  must  be  estimated  for  the  particular  application. 

If  for  a  particular  5,  Sm  does  not  satisfy  any  of  these  bounds,  then  convergence  is  not  guaran¬ 
teed  fay  the  above  theorems.  However,  the  following  general  technique  can  be  systematically  used 
to  overcome  the  problem,  provided  ||Sa||  <  1  in  some  norm.  Define  a  modified  iteration  operator 
S  by 

Mum 

S( u, t)  *  S(S  -  -  S(S(u,  I), t),  •  •  • ,  «),0- 

In  other  words,  5  is  obtained  by  iterating  5  k  times  with  t  fixed.  It  follows  that 


Therefore,  we  have 


>(S)  £  IIS.II  s  Iis.‘ll  £  IIS.II*. 

If  ||5«||  <  1,  then  a  large  enough  value  of  k  can  always  be  chosen  so  that  p(S)  or  ||5V||  satisfies 
one  of  the  bounds  in  Theorems  3.3  and  3.2.  For  efficiency  reasons,  k  should  be  chosen  to  be  the 
smallest  integer  such  that  the  largest  applicable  bound  is  satisfied. 

4.  Arc  length  Continuation 

In  this  section,  we  apply  the  results  of  the  last  section  to  an  important  application  area.  In 
arclength  continuation  methods  [2,  9,  12,  14],  G  represents  a  system  of  parameterized  nonlinear 
equations,  with  u  playing  the  role  of  the  main  variable,  t  the  parameters  and  N  represents  certain 
auxilliary  conditions.  We  shall  restrict  our  attention  to  path-following  continuation  where  m  =  1, 
although  the  algorithm  and  theory  developed  in  Sections  2  and  3  apply  to  other  related  problems 
as  well  (e.g.  augmented  systems  defining  singular  points  [1,  3,  11,  13,  15]). 

The  function  AT  is  usually  defined  in  terms  of  the  unit  tangent  (u,  t)  at  a  solution  (u,  t)  which 
is  the  solution  of: 

+  Gft  =  0 

IWII +  <’■=!■ 

We  shall  concentrate  on  two  typical  AT* s  that  are  widely  used  in  the  literature: 

AT1  =  u%(u  -  u0)  +  t0(t  -  to)  ~  6  a, 

N*=eJ (“I  1  <J  + 

where  (uo,  to)  Is  *  known  solution  on  the  solution  curve,  6a  is  a  continuation  step  and  t;-  is  the  j-th 
unit  vector.  For  more  details  the  reader  is  referred  to  [12]  for  Nl  and  [14]  for  AT*. 

We  shall  first  estimate  p{P)  and  ||P||  for  these  two  AT’s  and  then  apply  the  results  of  the  last 
section.  In  particular,  we  would  like  to  determine  the  conditions  under  which  Algorithm  ANM 
converges  if  S  is  convergent  for  G,  i.e  if  p(Sm)  <  1  . 

First,  we  need  the  following  elementary  result. 

Lemma  4.1. 

Proof.  Since  m  =  1,  P  is  a  rank  one  perturbation  of  /  and  thus  P  has  n  —  1  eigenvalues  equal  to 
1  and  one  eigenvalue  equal  to  1  +  j v,-1^ 

I 

For  Nl,  we  have  Nt  =  tij  and  Nt  =  tQ.  As  discussed  before,  if  5S  is  sufficiently  contractive, 
then  f 

u 

vm  -u,  «--r. 

It  follows  that  .  . 

,(#>)» 

If  tot  has  the  same  sign  as  u|u,  which  would  be  the  case  if  («o,to)  end  (u,t)  are  not  on  opposite 
sides  of  •  lrning  point,  then  we  have  p(P)  <  1.  If  tot  and  UqU  have  opposite  signs,  then  p[P)  >  1. 


*  imnptioa  tbit  Om  k  oontingnhr  nnrti  that  I  ft  0. 


But  note  that  the  term  tot  +  UqU  is  the  cosine  of  the  angle  #  between  the  unit  tangents  at  (u,  t) 
and  at  (uo,i0),  which  is  usually  kept  appreciably  above  zero  by  the  continuation  method.  Since 
|*o |  <  1  and  |t|  <  1,  we  have  in  this  case  p(P)  <  ~f.  In  particular,  as  6a  -*  0,  i  -*  0  and  hence 
p(P)  <  1.  Moreover,  if  5  is  sufficiently  contractive,  then  t>  tends  to  a  scalar  multiple  of  Nu  and 
this  implies  that  P  tends  to  being  a  normal  matrix.  As  for  ||P||,  we  have 

p=/-^V 

*0*  +  »0 O 


and  therefore 

l|P|la  -  1  +  P  =  2,00 ' 

COS  V 

As  6  s  — »  0,  ||P||F  <  2  for  p  —  2,oo.  If  P  is  normal,  we  have  the  tighter  bound  ||P||2  =  p(P)  <  1. 
Combining  the  above  estimates  of  p(P)  and  ||P||  with  the  results  of  Section  3,  we  obtain  the 
following  sufficient  conditions  for  the  convergence  of  Algorithm  ANM: 

Theorem  4.1.  For  N*,  as  6s  — ►  0,  Algorithm  ANM  converges  locally  if  any  one  of  the  following 
conditions  holds: 

I-  l|£«llr  <  |,  for  p  =  2  or  oo. 

2.  ||5a||2  <  1,  if  P  is  normal. 

3.  p(Sn)  <  |,  if  Sm  is  normal. 

4.  p(S%)  <  1,  if  P  and  S%  are  either  both  normal  or  simultaneously  diagonalizable. 

For  N2  we  have  (Nu,  Nt)  =  ej.  If  j  <  n  then  iVj  =  0  and  we  have  p(P)  =  1.  If  j  =  n  + 1  then 
AT«  *  0  and  we  also  have  p(P)  =  1.  Therefore  in  any  case,  p(P)  =  1.  However,  P  is  not  normal 
unless  v  is  a  multiple  of  ej.  Next  we  shall  estimate  ||P||.  First  note  that  if  j  =  n  +  1  then  P  —  I 
and  hence  ||P||  =  1  in  any  vector  induced  norm.  If  1  <  j  <  n,  then 


P  = 


where  {v)j  denotes  the  j-th  component  of  the  vector  v.  In  practice,  the  index  j  is  usually  chosen 
so  that  |(v)j-|  =  maxi<,<„  |(v},|  and  hence  ||P||oo  <  2  and  ||P||2  <  1  +  y/ri.  Combining  these 
estimates  of  p(P)  and  ||P||  with  the  results  of  Section  3  gives  the  following  sufficient  conditions  for 
the  convergence  of  Algorithm  ANM: 

Thaoram  4.S.  For  N2,  assuming  that  the  index  j  is  chosen  such  that  |(v),-|  =  maxi<,<„  |(«),-|, 
Algorithm  ANM  converges  if  any  one  of  the  following  conditions  holds: 

I’  Il^afloo  <  }■ 

lis-lb  <  r&- 

3.  ||5.||j  <  1,  if  P  is  normal. 

4.  f(Sm)  <  if  Su  is  normal. 

5.  p(Sm)  <  1,  if  P  and  SM  are  either  both  normal  or  simultaneously  diagonalizable. 


Many  of  the  conditions  in  Theorems  4.1  and  4.2  are  very  conservative.  For  AT1,  if  5a  is  normal 
and  reasonably  contractive,  then  it  is  most  likely  that  the  condition  p(Sv)  <  1  is  sufficient  because 
P  should  be  close  to  being  normal.  For  N2,  the  estimates  for  ||P||  are  especially  conservative  if  v 
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has  a  particularly  large  component.  In  fact,  if  is  equal  to  ey,  then  P  is  normal  and  ||P||p  =  I 
for  p  =  2  or  oo,  and  the  bounds  in  Theorem  4.2  all  become  1.  In  particular,  under  these  conditions, 
we  have  that,  for  both  Nl  and  N*,  the  condition  #(5«)  <  1  is  sufficient  for  the  convergence  of 
Algorithm  ANM  provided  SM  is  normal.  This  is  a  very  satisfactory  result  because  it  means  that  in 
practice  Algorithm  ANM  converges  if  S  is  convergent  for  G.  Therefore  it  can  be  applied  reliably 
to  a  large  class  of  problems  with  most  solvers  S  for  these  continuation  methods,  especially  in 
conjunction  with  the  technique  for  constructing  S.  The  algorithm  has  been  successfully  applied  to 
a  limited  number  of  small  continuation  problems. 

S.  Concluding  Remarks 

In  this  paper,  we  have  proposed  a  general  algorithm  for  solving  a  genera)  class  of  coupled 
nonlinear  systems.  It  is  especially  suitable  for  problems  for  which  efficient  solvers  exist  for  part  of 
the  system  but  not  for  the  whole.  The  algorithm  can  be  applied  in  a  modular  fashion  with  calls  to 
these  solvers  and  fully  exploits  the  efficiency  built  into  them.  The  local  convergence  analysis  shows 
that  if  the  solvers  are  sufficiently  contractive,  then  the  algorithm  converges  locally.  A  general 
technique  enables  one  to  construct  a  modified  S  that  ensures  convergence.  For  two  important 
arelength  continuation  methods,  we  show  that  under  mild  conditions,  the  algorithm  converges  if 
the  solver  is  convergent. 

It  is  known  that  the  Block  Elimination  Algorithm  may  be  unstable  near  a  solution  where  Gm  is 
singular  [5,  4].  If  C?^1  is  used  explicitly  in  S,  then  instability  is  to  be  expected  for  Algorithm  ANM 
as  well.  However,  it  may  be  possible  to  adapt  deflation  techniques  developed  in  [5,  6]  for  the  Block 
Elimination  Algorithm  to  Algorithm  ANM.  In  any  case,  for  problems  in  which  5  is  well-behaved, 
Algorithm  ANM  should  not  encounter  any  stability  problem. 
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